%--------------------------------------------------------------------------
% main.m
%--------------------------------------------------------------------------
% Original program by EJL (2012), adaptation by J. Mountjoy (2014)
%--------------------------------------------------------------------------
% Solves the optimal down and repayment problem for simulated borrowers for 
% a given parameter vector, computes moments, and draws scatter plots
%--------------------------------------------------------------------------


%--------------------------------------------------------------------------
% variable declarations
%--------------------------------------------------------------------------
clear; clc;

% cd(<folder directory>);

%parameters
global B T G Nsim mindwn ptilda dscale pscale deltaD deltaP APRUSE TERM;
global rho psi cval delta beta betap rra theta;
global saleeps downeps epsilon up ud;
global ystep ymax ymin ystates ynum y;
global lstep lmax lmin lstates lnum l;
global vstep vmax vmin vstates vnum v;

%--------------------------------------------------------------------------
% define discretized state space
%--------------------------------------------------------------------------
ymin=0.050; ystep=0.050; ymax=10.00; 
ystates=(ymin:ystep:ymax); ynum=size(ystates,2); %income states

lmin=0.100; lstep=0.100; lmax=2.000; 
lstates=(lmin:lstep:lmax); lnum=size(lstates,2); %monthly payment states (should be m in EJL notation)

vmin=-3.00; vstep=0.600; vmax=6.000; 
vstates=(vmin:vstep:vmax); vnum=size(vstates,2); %car valuation states

%--------------------------------------------------------------------------
% define parameters
%--------------------------------------------------------------------------
T = 40; G = T; B = 20; 
dscale=5; pscale=20; 
deltaD=0.100*dscale;  %Counterfactual for minimum down payment
deltaP=0.025*pscale;  %Counterfactual for car price

% offer terms from the sample
mindwn=1.000;       % sample average (rounded)
ptilda=11.000;      % sample average (rounded)
APRUSE=29.9;        % sample mode
TERM=42.0;          % sample mode

% imposed parameter values
delta = 0.880^((TERM/T)/12); 
beta  = 0.750^((TERM/T)/12); 
betap = beta; 
theta = 1.000; 
rra   = 0.999; 
rho   = 1.000; 
cval  = 0.000;
psi   = 0.000;

% calibrated parameter values (calibrated from proprietary data)
vavg = 2.9309;
vvar = 1.5348;
davg = 0.5012;
dvar = 2.0616;
corr =-0.0247;
sig2 = 0.2954;
phi  = 1.0295;

%--------------------------------------------------------------------------
% consumption utilities (YxLxV)
%--------------------------------------------------------------------------
y = repmat(ystates',[1,lnum,vnum]);
l = repmat(lstates ,[ynum,1,vnum]);
v = ones(1,1,vnum); v(1,1,:)=vstates; v=repmat(v,[ynum,lnum,1]);
up = CalcUtility(y,l);  % flow utility from having income y and car payment l, i.e. net money (y-l)
ud = CalcUtility(y,0);

%--------------------------------------------------------------------------
% simulate unobservables
%--------------------------------------------------------------------------
Nsim = 1*10^5;  % Number of simulated borrowers
%Nsim = 1*10^2;
stream = RandStream('mt19937ar','Seed',1); RandStream.setGlobalStream(stream); saleeps = randn(Nsim,1); %outside option
stream = RandStream('mt19937ar','Seed',2); RandStream.setGlobalStream(stream); downeps = randn(Nsim,1); %initial liquidity
stream = RandStream('mt19937ar','Seed',3); RandStream.setGlobalStream(stream); epsilon = randn(Nsim,T); %liquidity shocks

%--------------------------------------------------------------------------
% solve borrower's optimal purchase and repayment problem
%--------------------------------------------------------------------------

% simulate outside option (vout) and liquidity (y0,...,yT) variables
% v0: initial car valuation
% vi: initial car valuation grid index (# of steps from bottom of grid)
% y0: initial liquidity (income)
% yi: initial liquidity grid index
% yt: realized income each period (outcome of random walk)
% yindex: grid index for realized income each period
[v0 vi y0 yi yindex yt] = simulate_unobservables(vavg,vvar,davg,dvar,corr,sig2);

% solve repayment problem for all states (note: independent of d,p)
% P: income state transition matrix
% EVP: expected value of making payment
% EVD: expected value of defaulting
[P EVP EVD D] = optimal_repayment(sig2,delta,beta,psi,phi,theta);

% calculate optimal outcomes for each borrower, given v0 and y0
% dopt: optimal down payment
% simsale: indicator for purchase
% simxtra: additional down payment made in excess of minimum requirement
% simdef: indicator for default 
% simfrac: fraction of payments made

% first with standard $1000 minimum down and $11,000 car price
[dopt  lopt  uout ubuy  uopt  simsale  simless  simxtra  simdef simfrac] = optimal_down(mindwn,ptilda,v0,vi,y0,yi,yindex,EVP,EVD,D,cval,betap,theta,beta,phi);

%counterfactual: increase min down by deltaD = $500 and reduce price by
%deltaP = $500
[dopt2 lopt2 uout ubuy2 uopt2 simsale2 simless2 simxtra2 simdef2 simfrac2] = optimal_down(mindwn+deltaD,ptilda+deltaP,v0,vi,y0,yi,yindex,EVP,EVD,D,cval,betap,theta,beta,phi);

%--------------------------------------------------------------------------
% Summary statistics
%--------------------------------------------------------------------------

fprintf('Baseline prob of sale = %f\n', mean(simsale)); 
fprintf('Countefactual prob of sale = %f\n', mean(simsale2));

fprintf('Baseline prob of default = %f\n', mean(simdef(simsale == 1)));
fprintf('Counterfactual prob of default = %f\n', mean(simdef2(simsale == 1)));

fprintf('Baseline frac of payments = %f\n', mean(simfrac(simsale == 1)));
fprintf('Counterfactual frac of payments = %f\n\n\n', mean(simfrac2(simsale == 1)));

%--------------------------------------------------------------------------
% AC, MC_hat, and MS 
%--------------------------------------------------------------------------

%calculate average costs using EJL equation (24)

%present value of loan payments
%S: months of payments made
%D: downpayment (thousands)
%kappa: internal discount rate of 10% (EJL pg 1417)
%z: interest rate of 29.9%

S   = TERM*simfrac(simsale == 1);
D   = mindwn+simxtra(simsale == 1);
kappa = .1^(1/12);
z = (APRUSE/100)^(1/12);
numerator   = (1/kappa)*(1-exp(-kappa*S));
denominator = (1/z)*(1-exp(-z*TERM));
pv_payments = (numerator/denominator).*(ptilda-D);

S2  = TERM*simfrac(simsale2 == 1);
D2  = mindwn+deltaD+simxtra2(simsale2 == 1);
%kappa, z are the same
numerator2 = (1/kappa)*(1-exp(-kappa*S2));
%denominator is the same;
pv_payments2 = (numerator2/denominator).*(ptilda+deltaP-D2);

%present value of expected recovery prob of recovery X PV conditional on recovery from EJL Table 1
pv_recovery     = simdef(simsale == 1)*.78*1.597;
pv_recovery2    = simdef2(simsale2 == 1)*.78*1.597;

%company's cost of sale
%indirect_cost: indirect cost from EJL Table 5
%direct_cost: calculated suing ratio of direct cost to price from EJL Table 1
indirect_cost   = 2.481; 
direct_cost     = ptilda*(6.096/10.777);
cost_of_sale    = direct_cost + indirect_cost;

%q prime
q       = mean(simsale);
delta_q = mean(simsale2) - mean(simsale);
q_prime = delta_q/deltaP;

%average cost
ac_hat  = - mean(pv_payments + pv_recovery - cost_of_sale);

%marginal cost hat
ac2         = - mean(pv_payments2 + pv_recovery2 - cost_of_sale);
delta_ac    = ac2 - ac_hat;
ac_prime    = delta_ac/delta_q;
mc_hat      = q*(ac_prime) + ac_hat;

fprintf('q_i = %f\n',q);
fprintf('q_i_prime = %f\n', q_prime);
fprintf('ac_hat = %f\n', ac_hat);
fprintf('ac_prime = %f\n', ac_prime);
fprintf('mc_hat = %f\n', mc_hat);

%--------------------------------------------------------------------------
% plot 
%--------------------------------------------------------------------------

theta_grid = transpose(.05:.05:1);

for i = 1:length(theta_grid),
    
    t = theta_grid(i);
    
    ac(i) = ac_hat;
    mc(i) = (mc_hat - (1-t)*ac_hat)/t;
    ms(i) = -(q/q_prime)/t;
    
    social_markup(i) = t*ms(i) + (1-t)*(ac(i)-mc(i));
        
end

index = find(theta_grid == .2);
social_markup_2 = social_markup(index);

hold on;
plot(theta_grid,social_markup*1000,'-k','linewidth',1.5 );
ylabel('Social Markup (P-MC)','FontSize',12); xlabel('Market Power (\theta)','FontSize',12);
plot([0 1],[social_markup_2*1000 social_markup_2*1000],'--k',...
    [.2 .2], ylim,'--k','linewidth',.5 );
box on;

f_out = '../Figures/Social_markup.pdf';
hgsave(f_out)
print('-painters','-dpdf','-r300',f_out);


%% the rest of the codes are for outputs in EJL
%--------------------------------------------------------------------------
% Generate histograms for Figure A1
%--------------------------------------------------------------------------

% %output data for down payment histogram
% s = simsale;
% EDGES = [-Inf 0.00001:0.2:1.80001 Inf];
% modhist = histc(simxtra(s),EDGES)/sum(s);
% downout = [modhist(1:end-1); 1-mean(s)]
% 
% %output data for default time histogram
% simpmts = simfrac*T;
% EDGES = -0.5:2.0:41.5;
% modhist = histc(simpmts,EDGES)/sum(s);
% fracout = modhist

%--------------------------------------------------------------------------
% Generate scatter plot: baseline, prior to comparative statics
%--------------------------------------------------------------------------
%{
%create curves
yhi=6; use = y0>-2 & y0<yhi & v0<6 & v0>-2; use=y0>-99999; % flag outliers for exclusion
a = abs(ubuy-uout)<0.2;  % indicate the consumers on the margin of buying (first curve)
b = dopt>mindwn & abs(dopt-mindwn)<0.005; % indicate optimal down payment very close to minimum (second curve)
temp = [y0(a) v0(a) use(a)]; temp = sortrows(temp,1); % intial liq. and car value for these marginal types
yforcurve1 = temp(:,1); % liq. values along marginal line
vthreshold = temp(:,2); % car values along marginal line
use1 = temp(:,3)==1; % indicate consumers with non-outlier values and marginality
% Now do the same for the second curve, i.e. those on the margin of making
% natural down payment close to minimum value
temp = [y0(b) v0(b) use(b)]; temp = sortrows(temp,1);
yforcurve2 = temp(:,1); xthreshold = temp(:,2); use2 = temp(:,3)==1;

%scatter plots
%{
hold on;
scatter(y0(use & ~simsale),v0(use & ~simsale),2)  % scatter consumers who don't buy 
hold on;
scatter(y0(use & simless),v0(use & simless),2) % consumers who buy with min down
hold on;
scatter(y0(use & simsale & ~simless),v0(use & simsale & ~simless),2) % more than min down
%}

%threshold curves
hold on;
line(yforcurve1(use1),vthreshold(use1))
hold on;
line(yforcurve2(use2),xthreshold(use2))

%output data for scatter plot
%{
points = [y0 v0 simsale simless];
fid = fopen('scatter.txt','wt');
for i=1:Nsim,
    fprintf(fid,'%6.3f  ',points(i,:));
    fprintf(fid,'\n');
end;
fclose(fid);
%}

%--------------------------------------------------------------------------
% Generate scatter plot: comparative static one: increase mindown by $500
%--------------------------------------------------------------------------
{
%create curves
yhi=6; use = y0>-2 & y0<yhi & v0<6 & v0>-2; use=y0>-99999; % flag outliers for exclusion
a = abs(ubuy2-uout)<0.2;  % indicate the consumers on the margin of buying (first curve)
b = dopt2>mindwn+deltaD & abs(dopt2-(mindwn+deltaD))<0.005; % indicate optimal down payment very close to minimum (second curve)
temp = [y0(a) v0(a) use(a)]; temp = sortrows(temp,1); % intial liq. and car value for these marginal types
yforcurve1 = temp(:,1); % liq. values along marginal line
vthreshold = temp(:,2); % car values along marginal line
use1 = temp(:,3)==1; % indicate consumers with non-outlier values and marginality
% Now do the same for the second curve, i.e. those on the margin of making
% natural down payment close to minimum value
temp = [y0(b) v0(b) use(b)]; temp = sortrows(temp,1);
yforcurve2 = temp(:,1); xthreshold = temp(:,2); use2 = temp(:,3)==1;

%scatter plots
figure;
hold on;
scatter(y0(use & ~simsale2),v0(use & ~simsale2),2)  % scatter consumers who don't buy 
hold on;
scatter(y0(use & simless2),v0(use & simless2),2) % consumers who buy with min down
hold on;
scatter(y0(use & simsale2 & ~simless2),v0(use & simsale2 & ~simless2),2) % more than min down

%threshold curves
hold on;
line(yforcurve1(use1),vthreshold(use1))
hold on;
line(yforcurve2(use2),xthreshold(use2))
%}

%--------------------------------------------------------------------------
% Generate scatter plot: comparative static two: increase price by $2000
%--------------------------------------------------------------------------

% %create curves
% yhi=6; use = y0>-2 & y0<yhi & v0<6 & v0>-2; use=y0>-99999; % flag outliers for exclusion
% a = abs(ubuy3-uout)<0.2;  % indicate the consumers on the margin of buying (first curve)
% b = dopt3>mindwn & abs(dopt3-mindwn)<0.005; % indicate optimal down payment very close to minimum (second curve)
% temp = [y0(a) v0(a) use(a)]; temp = sortrows(temp,1); % intial liq. and car value for these marginal types
% yforcurve1 = temp(:,1); % liq. values along marginal line
% vthreshold = temp(:,2); % car values along marginal line
% use1 = temp(:,3)==1; % indicate consumers with non-outlier values and marginality
% % Now do the same for the second curve, i.e. those on the margin of making
% % natural down payment close to minimum value
% temp = [y0(b) v0(b) use(b)]; temp = sortrows(temp,1);
% yforcurve2 = temp(:,1); xthreshold = temp(:,2); use2 = temp(:,3)==1;
% 
% %scatter plots
% figure;
% hold on;
% scatter(y0(use & ~simsale3),v0(use & ~simsale3),2)  % scatter consumers who don't buy 
% hold on;
% scatter(y0(use & simless3),v0(use & simless3),2) % consumers who buy with min down
% hold on;
% scatter(y0(use & simsale3 & ~simless3),v0(use & simsale3 & ~simless3),2) % more than min down
% 
% %threshold curves
% hold on;
% line(yforcurve1(use1),vthreshold(use1))
% hold on;
% line(yforcurve2(use2),xthreshold(use2))


%--------------------------------------------------------------------------
% Figure A3(a): purchase utility and optimal down versus price
%--------------------------------------------------------------------------
%{
numpts = 21;
values = zeros(numpts,5);
simulated_apps = [];
simulated_sales = [];
for i=1:numpts,
    
    %calculate simulated outcomes conditional on sale
    [dopt lopt uout ubuy uopt simsale simless simxtra simdef simfrac] = optimal_down(mindwn,i,v0,vi,y0,yi,yindex,EVP,EVD,D,cval,betap,theta,beta,phi);
    x1 = mean(simsale);
    x2 = mean(simless(simsale));
    x3 = mean(dopt(simsale));
    x4 = mean(simfrac);
    x5 = mean(simdef);
    values(i,:)=[x1 x2 x3 x4 x5];
    %{
    %caluclate simulated outcomes not conditional on sale
    temp_d = -99; %no minimum down restriction
    [dopt lopt uout ubuy uopt simsale simless simxtra simdef simfrac] = optimal_down(temp_d,i,v0,vi,y0,yi,yindex,EVP,EVD,D,cval,betap,theta,beta,phi);
    temp = dopt; temp(temp<0)=0;
    x6 = mean(temp);
    x7 = mean(simfrac);
    x8 = mean(simfrac);
    values(i,:)=[x1 x2 x3 x4 x5 x6 x7 x8];
    %}
    %create data for reduced-form regressions
    new_apps = [y0 i*ones(Nsim,1) mindwn*ones(Nsim,1) simsale simless];
    new_sales = [new_apps(simsale,:) dopt(simsale) simdef simfrac];
    simulated_apps = [simulated_apps; new_apps];
    simulated_sales = [simulated_sales; new_sales];
    
end;
asdf
%--------------------------------------------------------------------------
% Output data
%--------------------------------------------------------------------------
fid = fopen('simoutcomes_vs_price.txt','wt');
for i=1:numpts,
    fprintf(fid,'%10.6f ',values(i,:));
    fprintf(fid,'\n');
end;
fclose(fid);
fid = fopen('simulated_apps_p.txt','wt');
for i=1:size(simulated_apps,1),
    fprintf(fid,'%10.6f,',simulated_apps(i,:));
    fprintf(fid,'\n');
end;
fclose(fid);
fid = fopen('simulated_sales_p.txt','wt');
for i=1:size(simulated_sales,1),
    fprintf(fid,'%10.6f,',simulated_sales(i,:));
    fprintf(fid,'\n');
end;
fclose(fid);
asdf
%}
%--------------------------------------------------------------------------
% Figure A3(b): purchase utility and optimal down versus min down
%--------------------------------------------------------------------------
%{
numpts = 1;
values = zeros(numpts,5);
simulated_apps = [];
simulated_sales = [];
for i=1:numpts,

    %calculate simulated outcomes conditional on sale
    [dopt lopt uout ubuy uopt simsale simless simxtra simdef simfrac] = optimal_down(0.1*(i-1),ptilda,v0,vi,y0,yi,yindex,EVP,EVD,D,cval,betap,theta,beta,phi);
    x1 = mean(simsale);
    x2 = mean(simless(simsale));
    x3 = mean(dopt(simsale));
    x4 = mean(simfrac);
    x5 = mean(simdef);
    values(i,:)=[x1 x2 x3 x4 x5];
    %{
    %caluclate simulated outcomes not conditional on sale
    temp_d = -99; %no minimum down restriction
    [dopt lopt uout ubuy uopt simsale simless simxtra simdef simfrac] = optimal_down(temp_d,ptilda,v0,vi,y0,yi,yindex,EVP,EVD,D,cval,betap,theta,beta,phi);
    temp = dopt; temp(temp<0)=0;
    x6 = mean(temp);
    x7 = mean(simfrac);
    x8 = mean(simfrac);
    values(i,:)=[x1 x2 x3 x4 x5 x6 x7 x8];
    %}
    %create data for reduced-form regressions
    new_apps = [y0 ptilda*ones(Nsim,1) 0.1*(i-1)*ones(Nsim,1) simsale simless];
    new_sales = [new_apps(simsale,:) dopt(simsale) simdef simfrac];
    simulated_apps = [simulated_apps; new_apps];
    simulated_sales = [simulated_sales; new_sales];
    
end;

%--------------------------------------------------------------------------
% Output data
%--------------------------------------------------------------------------
fid = fopen('simoutcomes_vs_mindp.txt','wt');
for i=1:numpts,
    fprintf(fid,'%10.6f ',values(i,:));
    fprintf(fid,'\n');
end;
fclose(fid);
fid = fopen('simulated_apps_d.txt','wt');
for i=1:size(simulated_apps,1),
    fprintf(fid,'%10.6f,',simulated_apps(i,:));
    fprintf(fid,'\n');
end;
fclose(fid);
fid = fopen('simulated_sales_d.txt','wt');
for i=1:size(simulated_sales,1),
    fprintf(fid,'%10.6f,',simulated_sales(i,:));
    fprintf(fid,'\n');
end;
fclose(fid);
%}
%--------------------------------------------------------------------------
% end of program
%--------------------------------------------------------------------------
